home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / histogram / pdf2d.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  4.3 KB  |  187 lines
  1. /* histogram/pdf2d.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_errno.h>
  23. #include <gsl/gsl_histogram.h>
  24. #include <gsl/gsl_histogram2d.h>
  25.  
  26. #include "find.c"
  27.  
  28. int
  29. gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * p,
  30.                 double r1, double r2,
  31.                 double *x, double *y)
  32. {
  33.   size_t k;
  34.   int status;
  35.  
  36. /* Wrap the exclusive top of the bin down to the inclusive bottom of
  37.    the bin. Since this is a single point it should not affect the
  38.    distribution. */
  39.  
  40.   if (r2 == 1.0)
  41.     {
  42.       r2 = 0.0;
  43.     }
  44.   if (r1 == 1.0)
  45.     {
  46.       r1 = 0.0;
  47.     }
  48.  
  49.   status = find (p->nx * p->ny, p->sum, r1, &k);
  50.  
  51.   if (status)
  52.     {
  53.       GSL_ERROR ("cannot find r1 in cumulative pdf", GSL_EDOM);
  54.     }
  55.   else
  56.     {
  57.       size_t i = k / p->ny;
  58.       size_t j = k - (i * p->ny);
  59.       double delta = (r1 - p->sum[k]) / (p->sum[k + 1] - p->sum[k]);
  60.       *x = p->xrange[i] + delta * (p->xrange[i + 1] - p->xrange[i]);
  61.       *y = p->yrange[j] + r2 * (p->yrange[j + 1] - p->yrange[j]);
  62.       return GSL_SUCCESS;
  63.     }
  64. }
  65.  
  66. gsl_histogram2d_pdf *
  67. gsl_histogram2d_pdf_alloc (const size_t nx, const size_t ny)
  68. {
  69.   const size_t n = nx * ny;
  70.  
  71.   gsl_histogram2d_pdf *p;
  72.  
  73.   if (n == 0)
  74.     {
  75.       GSL_ERROR_VAL ("histogram2d pdf length n must be positive integer",
  76.             GSL_EDOM, 0);
  77.     }
  78.  
  79.   p = (gsl_histogram2d_pdf *) malloc (sizeof (gsl_histogram2d_pdf));
  80.  
  81.   if (p == 0)
  82.     {
  83.       GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf struct",
  84.             GSL_ENOMEM, 0);
  85.     }
  86.  
  87.   p->xrange = (double *) malloc ((nx + 1) * sizeof (double));
  88.  
  89.   if (p->xrange == 0)
  90.     {
  91.       free (p);        /* exception in constructor, avoid memory leak */
  92.  
  93.       GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf xranges",
  94.             GSL_ENOMEM, 0);
  95.     }
  96.  
  97.   p->yrange = (double *) malloc ((ny + 1) * sizeof (double));
  98.  
  99.   if (p->yrange == 0)
  100.     {
  101.       free (p->xrange);
  102.       free (p);        /* exception in constructor, avoid memory leak */
  103.  
  104.       GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf yranges",
  105.             GSL_ENOMEM, 0);
  106.     }
  107.  
  108.   p->sum = (double *) malloc ((n + 1) * sizeof (double));
  109.  
  110.   if (p->sum == 0)
  111.     {
  112.       free (p->yrange);
  113.       free (p->xrange);
  114.       free (p);        /* exception in constructor, avoid memory leak */
  115.  
  116.       GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf sums",
  117.             GSL_ENOMEM, 0);
  118.     }
  119.  
  120.   p->nx = nx;
  121.   p->ny = ny;
  122.  
  123.   return p;
  124. }
  125.  
  126. int
  127. gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * p, const gsl_histogram2d * h)
  128. {
  129.   size_t i;
  130.   const size_t nx = p->nx;
  131.   const size_t ny = p->ny;
  132.   const size_t n = nx * ny;
  133.  
  134.   if (nx != h->nx || ny != h->ny)
  135.     {
  136.       GSL_ERROR ("histogram2d size must match pdf size", GSL_EDOM);
  137.     }
  138.  
  139.   for (i = 0; i < n; i++)
  140.     {
  141.       if (h->bin[i] < 0)
  142.     {
  143.       GSL_ERROR ("histogram bins must be non-negative to compute"
  144.              "a probability distribution", GSL_EDOM);
  145.     }
  146.     }
  147.  
  148.   for (i = 0; i < nx + 1; i++)
  149.     {
  150.       p->xrange[i] = h->xrange[i];
  151.     }
  152.  
  153.   for (i = 0; i < ny + 1; i++)
  154.     {
  155.       p->yrange[i] = h->yrange[i];
  156.     }
  157.  
  158.   {
  159.     double mean = 0, sum = 0;
  160.  
  161.     for (i = 0; i < n; i++)
  162.       {
  163.     mean += (h->bin[i] - mean) / ((double) (i + 1));
  164.       }
  165.  
  166.     p->sum[0] = 0;
  167.  
  168.     for (i = 0; i < n; i++)
  169.       {
  170.     sum += (h->bin[i] / mean) / n;
  171.     p->sum[i + 1] = sum;
  172.       }
  173.   }
  174.  
  175.   return GSL_SUCCESS;
  176. }
  177.  
  178.  
  179. void
  180. gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * p)
  181. {
  182.   free (p->xrange);
  183.   free (p->yrange);
  184.   free (p->sum);
  185.   free (p);
  186. }
  187.